import numpy as np
import pandas as pd
from sympy import *

data_2 = pd.read_csv('data2.csv')
data_3=pd.read_csv('data3.csv')
time = np.array(data_2['time'])
# time = np.array(data_3['time'])
ti = np.array([int(t[:2]) * 60 + int(t[3:]) for t in time])
x = np.array(data_2['x'])
y = np.array(data_2['y'])

omega=np.array([atan(yi/xi) for (xi,yi) in zip(x,y)])
N0=79.6764+0.2422*(2015-1985)-int((2015-1985)/4)
x0 = np.array([sqrt(xi * xi + yi * yi) for (xi, yi) in zip(x, y)])
c2 = 116

# 通过误差的约束 搜索枚举出可能的地点和日期
res=[]
for hi in np.arange(1, 4, 0.5):
    theta = np.array([atan(hi / xi) for xi in x0])
    for Ni in np.arange(1,366,1):
        t=np.array([asin((1.686*10**12*sqrt(6.315-cos(0.0172*(Ni-173))**2)*cos(th)*sin(o))/(-4.24*10**12+6.71*10**11*cos(0.0172*(Ni-173))**2)) for (o,th) in zip(omega,theta)])
        beta=asin(0.398*cos(0.9856*(Ni-173)*3.1415/180))
        t0=np.array([((ti*180)/(15*3.14)+12)*60 for ti in t])
        c1=np.array([c2-(tii-t0i)/4 for (tii,t0i) in zip(ti,t0)])

        for psi in np.arange(-90,90,0.5):
            wucha=np.array([sin(th)-sin(beta)*sin(psi)-cos(beta)*cos(psi)*cos(ti) for (th,ti) in zip(theta,t)])

            print(hi, Ni, psi)
            if abs(wucha.mean())<0.00001:
                print(float(wucha.mean()))
                print(hi,Ni,c1.mean(),psi)
                res.append([hi,Ni,c1.mean(),psi])

print(res)